from math import pi
from numpy import *
import matplotlib.pyplot as plt
from cplinalg import *
def rechighcut(x,wc,sr,flen):
	xlf=fft.rfft(x,flen)
	dw=2*pi/sr/flen
	wb=round(float(wc)/dw)
	if wb<xlf.size:
		xlf[wb:]=0
	y=fft.irfft(xlf,flen)
	return y
		
def recbndpass(x,lc,hc,sr,flen):
	xlf=fft.rfft(x,flen)
	dw=2*pi/sr/flen
	hb=round(float(hc)/dw)
	lb=round(float(lc)/dw)
	if hb<xlf.size:
		xlf[hb:]=0
	xlf[0:lb+1]=0
	y=fft.irfft(xlf,flen)
	return y

def hrft(x):
	N=x.size
	dx=(roll(x,-1)-roll(x,1))/2.0
	dx[0]=x[1]-x[0]
	dx[-1]=x[-1]-x[N-2]
	targ=vstack([x.reshape(N,1),dx.reshape(N,1)])
	t=arange(0,N)
	w=arange(0,2*N)
	ind=dot(t.reshape(N,1),w.reshape(1,2*N))
	Mu=exp(1j*pi/N*ind)
	wd=w
	wd[N+1:]=-w[N-1:0:-1]
	D=diag(wd*pi*1j/N,0)
#	print D
	Mb=dot(Mu,D)
	M=vstack([Mu,Mb])
#	sv=linalg.svd(M)[1]
#	print sv
	y=tiklstsq(M,targ,1e-3).reshape(2*N,)
#	y=linalg.lstsq(M,targ,1e-12)[0].reshape(2*N,)
	return y
	
	
		


if __name__=="__main__":
	L=256
	t=arange(0,L)
	x1=sin(5*pi/L*t)
	x2=sin(7*pi/L*t)
	x=x1+x2
	y=hrft(x)
	yc=roll(y,256)
	yy=fft.fft(x)
	scl=abs(yy[0])/abs(y[0])
	yyc=roll(yy,128)
	fy=arange(-128,128,0.5)
	fyy=arange(-128,128)
#	print fy.shape
#	print fyy.shape
#	print yyc.shape
#	fl=int(20*L)
#       wc=pi/L
#	print wc
#	y=rechighcut(x,wc,1,fl)
#	xl=zeros(y.size)
#	xl[0:x.size]=x
	plt.plot(fy,log(abs(yc)*scl))
#	plt.plot(fyy,log(abs(yyc)))
#	print abs(y)
#	plt.plot(imag(y))
	plt.show()
